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I outline the theory of relativistic charged-particle motion in the magnetosphere in a way suitable 
for undergraduate courses. I discuss particle and guiding center motion, derive the three adiabatic 
invariants associated with them, and present particle trajectories in a dipolar field. I provide twelve 
computational exercises that can be used as classroom assignments or for self-study. Two of the 
exercises, drift-shell bifurcation and Speiser orbits, are adapted from active magnetospheric research. 
The Python code provided in the supplement can be used to replicate the trajectories and can be 
easily extended for different field geometries. 



I. INTRODUCTION 

Ions and electrons trapped in the Earth's magnetic 
field may affect our technology and our daily lives in sig- 
nificant ways. Energetic plasma particles may penetrate 
satellites and disable them temporarily or permanently. 
They can also pose serious health hazards for astronauts 
in space. Spectacles like the aurora are created by parti- 
cles that enter the Earth's atmosphere at polar regions; 
on the other hand, aircraft personnel and frequent flyers 
may accumulate a significant dose of radiation due to the 
same particles. 1 All of these effects are enhanced at pe- 
riods of solar maximum, the next one being expected to 
happen between 2012 and 2014. Occasional extreme so- 
lar events may induce currents in the ionosphere, which 
in turn induce significant currents on power lines, causing 
power outages. 2 Such events can also disrupt communi- 
cations, radio and GPS. Thus, understanding and pre- 
dicting the processes in the Earth's magnetosphere have 
practical importance. 

This paper aims to outline one of these processes, 
charged-particle motion and associated adiabatic invari- 
ants, for physics students and instructors who wish to use 
it in lectures. The emphasis is on numerical computation 
and visualization of trajectories. For a more comprehen- 
sive discussion advanced texts on plasma physics 3 - - — can 
be consulted. 

Other authors have also suggested using topics from 
plasma research to enhance undergraduate curriculum. 
Lopez^ provides examples of how space physics can be in- 
corporated in undergraduate electromagnetism courses, 
and McGuire 8 shows how computer algebra systems can 
be used to follow particle trajectories in electric dipole 
and (separately) in magnetic dipole fields. Photographs 
of plasma experiments, such as those provided by Hug- 
gins and Lelek 9 and by the UCLA Plasma Lab web site^ 
are also helpful for understanding space plasma behavior. 

Figure [T] shows a schematic description of the Earth's 
magnetosphere, which is the region in space where the 
magnetic field of the Earth is dominant. Charged par- 
ticles trapped in the magnetosphere form the radiation 
belts, the plasmasphere, and current systems such as the 
ring current, tail current, and field-aligned currents. 

The Earth radius R e (6378.137 km) is a natural length 
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FIG. 1. A schematic view of the Earth's magnetosphere. The 
solar wind comes from the left. (Courtesy of Kenneth R. 
Lang, 11 reproduced with permission.) 



scale for the magnetosphere. Near the Earth, up to 3- 
4R e , the field can be very well approximated with the 
field of a dipole. However, at larger distances, the effects 
of the solar wind cause significant deviations from the 
dipole. 

The solar wind is a stream of plasma carrying mag- 
netic field from the Sun. When the solar wind encounters 
the Earth's magnetosphere, the two systems do not mix. 
This is because of the "frozen-in flux" condition 12 which 
dictates that plasma particles stay attached to magnetic 
field lines, except at special locations such as polar cusps. 
The solar wind influences the magnetosphere by apply- 
ing mechanical and magnetic pressure on it, compressing 
it earthward on the side facing the Sun (the "dayside"). 
This compression is stronger when the Sun is more ac- 
tive. On the opposite side (the "nightside" ) , the field is 
extended over a very large distance, forming the mag- 
netotail. Wolf 13 provides a review of the complex and 
time-dependent interactions between magnetic fields, in- 
duced electric fields and plasma populations. 

The Van Allen radiation belts form a particularly sig- 
nificant plasma population due to their high energy and 
their proximity to Earth. They can be found from 
1000km above the ground up to a distance of 6R e . These 
belts are composed of electrons with energies up to sev- 
eral MeVs and of protons with energies up to several 
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hundred Me Vs. The dynamics of these particles is the 
main focus of this paper. 

This paper is organized as follows: Section |TT] intro- 
duces the relativistic equation of motion for a particle 
in an electric and magnetic field and describes the cy- 
clotron, bounce and drift motions. It also shows some 
typical particle trajectories under the dipolar magnetic 
field, approximating the Earth's field. Section IIIII intro- 
duces the concept of adiabatic invariants and derives the 
first adiabatic invariant associated with the particle mo- 
tion. Section HVl gives the approximate equations of mo- 
tion for the guiding center of a particle, obtained by aver- 
aging out the cyclotron motion. Section [V] presents and 
derives the second and third invariants associated with 
the bounce and drift motions, respectively. Section |VI] 
lists some exercises building on the concepts described 
in the paper. Two of these exercises describe non-dipole 
fields that are used for modeling different regions of the 
magnetosphere. 



II. PARTICLE TRAJECTORY IN DIPOLAR 
MAGNETIC FIELD 

The motion of a particle with charge q and mass m in 
an electric field E and magnetic field B is described by 
the Newton- Lor entz equation: 



d(7mv) 
dt 



qE(r) +gvx B(r). 



(1) 



Here 7 = (1 — v 2 / 'c 2 ) - ^ 1 / 2 ) is the relativistic factor and v 
is the particle speed. 

Suppose that E = 0. Then, because of the cross prod- 
uct, the acceleration of the particle is perpendicular to 
the velocity at all times, so the speed of the particle (and 
the factor 7) remains constant. Further suppose that the 
magnetic field is uniform. Then, particles move on he- 
lices parallel to the field vector. The circular part of this 
motion is called the "cyclotron motion" or the "gyromo- 
tion". The "cyclotron frequency" Q and the "cyclotron 
radius" p are respectively given by 
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(2) 
(3) 



where B — |B| is the uniform field strength and v± is 
the component of the velocity perpendicular to the field 
vector. If there are not any other forces, the parallel com- 
ponent of the velocity remains constant. The combined 
motion traces a helixJ^ 

If the electric field is not zero, we can write it as 
E = E^ + En, where E|| is parallel to B, and E^ is 
perpendicular to it. If E|| 7^ 0, particles accelerate with 
qE\\/m along the field line and they are rapidly removed 
from the region. Therefore, the existence of a trapped 
plasma population implies that the parallel electric field 
must be negligible. 



The perpendicular component of the electric field will 
move particles with an overall drift velocity, known as 
the E-cross-B-drift, which is perpendicular to both field 
vectors: 



EiXB 
B 2 ' 



(4) 



The particle will move with with the velocity v#, plus the 
cyclotron motion described above. The drift velocity v# 
is independent of particle mass and charge. Therefore, 
in an inertial frame moving with v#, the E-cross-B-drift 
will vanish for all types of particles. 

For the remainder of this paper we take E = 0, that 
is, the acceleration due to the electric field is not taken 
into consideration. This is not because electric fields are 
unimportant; on the contrary, they play an important 
role in the complex dynamics of plasmas. The first reason 
for leaving out electric field effects is described above: 
If the field is uniform and constant, we can transform 
to another frame that cancels it. Even if the field is 
nonuniform and time-dependent, electric drifts can be 
vectorially added to magnetic drifts in order to obtain 
the overall drift. Drift velocities due to different fields 
are independent. 

The second reason is the need for simplicity; a static 
magnetic field provides sufficient real-life context for the 
discussion of guiding-center and adiabatic invariant con- 
cepts in general-purpose lectures. The final reason is that 
this paper focuses on the region occupied by radiation 
belts, and in this region the magnetic term of (pQ) is the 
dominant forced 

Now we consider motion under the influence of a mag- 
netic dipole. The field Bdi p (r) of a magnetic dipole with 
moment vector M at location r is given by>^ 



B dip (r) 
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[3(M • r)r - M] , 



(5) 



where r = xx + yy + zz, r = |r| and r = r/r. For 
Earth, we take M = — Mz, antiparallel to the z-axis, 
because the magnetic north pole is near the geographic 
south poleJ£ 

At the magnetic equator (x = lR ei y — z = 0) the field 
strength is measured to be Bo = 3.07 x 10 _5 T. Substi- 
tution shows that hqM/At: = BqR^. Then in Cartesian 
coordinates, the field is given by: 

B dip = [Sxzxl + iyzy + (2z 2 - x 2 - y 2 )z] . (6) 

Figure [2] shows trajectories of two protons with lOMeV 
kinetic energy, a typical energy for radiation belts. 13 
The trajectories are calculated with the SciPy module 
using the Python language^ One proton is started at 
(2i? e ,0,0) and the other at (4R e ,0,0). Both start with 
an equatorial pitch angle (angle between the velocity 
and field vectors) a eq = 30° so that v y o = vsina e q, 
v z o = vcoso; eq and v x o = 0. Both are followed for 120 
seconds. 
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the drift motion period t<± are approximately given as£ 
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FIG. 2. (Color online) Trajectories of two lOMeV protons in 
the Earth's dipole field. The dipole moment is in the — z di- 
rection. Both panels show the same trajectories from different 
viewing angles. 



The motion is again basically helical, but the nonuni- 
formity of the field introduces two additional modes of 
motion on large spatial and temporal scales. These are 
called "the bounce motion" and "the drift motion". 

The bounce motion proceeds along the field line that 
goes through the helix (the "guiding line"). The motion 
slows down as it moves toward locations with a stronger 
magnetic field, reflecting back at "mirror points". The 
bounce motion is much slower than the cyclotron motion. 

The drift motion takes the particle across field lines 
(perpendicular to the bounce motion). In general, drift 
motion is faster at larger distances, as observed in Fig- 
ure [21 Particles in dipole-like fields are trapped on closed 
"drift shells" as long as they are not disturbed by colli- 
sions or interactions with EM waves. The drift motion is 
much slower than the bounce motion. 

Under a dipolar field, the bounce motion period Tb and 
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where Rq is the equatorial distance to the guiding line 
and a eq is the equatorial pitch angle. Both approxima- 
tions have an error of about 0.5%. 



III. THE FIRST ADIABATIC INVARIANT 

If the parameters of an oscillating system are varied 
very slowly compared to the frequency of oscillations, 
the system possesses an "adiabatic invariant", a quan- 
tity that remains approximately constant. In the Hamil- 
tonian formalism, the adiabatic invariant is the same as 
the action variable)^ 



J 



p(E) dq, 



H=E 



(9) 



where g, p are canonical variables and the integral is eval- 
uated over one cycle of motion satisfying H(p : q,t) = E. 
The integral should be evaluated at "frozen time", that 
is, the slowly varying parameter is considered constant 
during the integration cycle. 

There are three separate periodic motions of a charged 
particle in a dipole-like magnetic field. This means there 
are three adiabatic invariants for the particle's motion. 

The canonical momentum for a charged particle in a 
magnetic field with vector potential A is 7m v + qA. To 
obtain the first adiabatic invariant Ji, we integrate the 
canonical momentum over one cycle of the cyclotron or- 
bit: 



J\ = j> (jmv + qA) • di. 



(10) 



Here d£ is the line element of the particle trajectory. 
Even though the path does not exactly close, we eval- 
uate the integral as if it does. Stokes' Theorem states 
that: 



A-dl 



V x A • dcr, 



(11) 



where the right-hand side integral is taken over the sur- 
face bounded by the closed path of the left-hand side 
integral. Using Stokes' theorem with B = V x A, the 
integral J\ takes the form: 



Ji = 27rjmpv± — qnp B 



9 99 
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qB 



(12) 
(13) 



The second equation follows from substituting p from (|3]) . 
It is assumed that the cycle is sufficiently small so that B 
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is considered uniform over the area bounded by the gyro- 
motion. This assumption is essential for the existence of 
adiabatic invariants. The area element da is antiparallel 
to B because of the sense of gyration of the particles; 
hence the negative sign of the second term. 
Customarily, one takes 



2 2 

2D 



(14) 



as the first invariant, which differs from J\ only in some 
constant factors. The parameter \i is called the "mag- 
netic moment" because it is equal to the magnetic mo- 
ment of the current generated by the particle moving on 
this circular path. 

The adiabatic invariance of \i explains the existence of 
mirror points. As the particle moves along the field line 
toward points with larger £?, the perpendicular speed v± 
must increase in order to keep \i constant. However, v± 
can be at most the constant total speed. The particle 
stops at the point where B = B m = j 2 mv 2 /2/j J and falls 
back (see Section HV]) . 

The left panel of Figure [3] shows the magnetic moment 
\i for the two protons shown in Figure [2j The values 
are oscillating with the local cyclotron frequency because 
instantaneous values of v± and B are used. The actual 
adiabatic invariant is the average of these oscillations and 
it is constant in time. 

The top right panel shows the oscillations of \i vs. time 
for a shorter interval. Comparison with the z vs. time 
plot below, it can be seen that these oscillations are cor- 
related with the bounce motion. Oscillations of \i have a 
small amplitude near the mirror point because there the 
parallel motion slows down and the overall motion be- 
comes more adiabatic. Near the equatorial plane (z = 0) 
parallel motion is fastest, the motion is less adiabatic, 
and /i oscillates with a larger amplitude. 

The proper way of calculating the first invariant would 
remove all oscillations: After following the full trajectory, 
find the times {U} where v x = (or v y = 0) by search- 
ing along discrete path points and by interpolating. The 
difference between any successive time points is half a 
gyroperiod. Then take the average of \i over that time 
interval using the path points. Repeating this procedure 
for all time intervals, we get a constant set of \i values 
(apart from numerical errors). 



IV. THE GUIDING-CENTER EQUATIONS 

The "guiding center" is the geometric center of the 
cyclotron motion. If the magnetic field is uniform the 
guiding center moves with constant velocity parallel to 
the field line. In nonuniform field geometries, there is a 
sideways drift in addition to the motion along the field 
lines, as seen in the dipole example in Figure [2j 

Calculation of the guiding center motion requires that 
the motion is helical in the smallest scale, and that the 



field does not change significantly within a cyclotron ra- 
dius. This condition can be expressed as 



B 



VB\ 



(15) 



The magnetic moment is an adiabatic invariant under 
this condition. 

Northrop 19 and Walt 5 give detailed derivations of the 
equations of guiding- center motion. In order to derive the 
acceleration of the guiding center, the particle position r 
is substituted with: 



r = R + p 



(16) 



where R is the position of the guiding center. The vector 
p lies on the plane perpendicular to the field, oscillates 
with the cyclotron frequency, and its length is equal to 
the cyclotron radius. 

Assuming that the cyclotron radius is much smaller 
than the length scale of the field, we can expand B(r) 
around R to first order in a Taylor series: 



B(r) «B(R) + (p. V)B 



(17) 



This expansion is substituted into the Newton-Lorenz 
equation (pQ) and the equation is averaged over a cycle, 
eliminating rapidly oscillating terms containing p and its 
derivatives. The resulting acceleration of the guiding cen- 
ter is given by: 
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(18) 



Taking the dot product of both sides of the equation 
with b, the local magnetic field direction, will yield the 
equation of motion along the field line. The first term 
becomes identically zero because it is perpendicular to 
the field vector. Then: 



dv 



qp 2 tt 



-li = -^lb-VB(K) 



dt 



27m 



(19) 



where v\\ is the speed along the field line. Replacing 
qp 2 Q/ (27m) = /i/(7 2 m) and defining s as the distance 
along the field line, this equation can be written as: 



(to,, 
dt 



d^s 
dt 2 



d fiB(s) 
ds 7 2 m 



(20) 



Here B(s) is the field strength along the field line. The 
factors \i and 7 can be taken inside the derivative be- 
cause they are constants. This expression shows that the 
quantity pB(s)/(j 2 m) acts like a potential energy in the 
parallel direction. The negative sign indicates that the 
parallel motion is accelerated toward regions with smaller 
field strength. 

The motion of the guiding center perpendicular to field 
lines can be determined by taking the cross product of 
Eq. ([15]) with the field direction vector. The resulting 
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FIG. 3. (Color online) Left: The instantaneous values of the magnetic moment for the particle orbits shown in Fig. [2] Lower 
curve is for the proton starting at 2 R e distance, upper curve for the proton starting at 4 R e . The white horizontal line shows 
the average value. Top right: A close-up view of the instantaneous magnetic moment up to time 6s. Bottom right: The 
z-coordinate up to time 6s. 



equation is then iterated to obtain an approximate solu- 
tion for the drift velocity across field lines. 



dRj 



7m 
2qB 2 
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)bxVB 



(21) 



This drift velocity is actually the sum of two separate 
drift velocities: The gradient drift that arises from the 
nonuniformity of the magnetic field, and the curvature 
drift that occurs because the field lines are curved. For 
an example of pure gradient drift motion, see Exercise [2] 
in Section EH 

Equation ([2T]) shows that electrons and ions drift in 
opposite directions. This creates a net current around 
the Earth, called "the ring current" . 

Gradient and curvature drifts are the only drifts seen 
in static magnetic fields. External electric fields, external 
forces such as gravity, and time dependent fields create 
additional drift velocities^ 

Combining these, we obtain the following equations of 
motion for the guiding center: 



dR 


■ymv 2 




2qB 2 
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fi 


dt 


7 2 m 




bxVB 



b-VB 



(22) 
(23) 



These equations are more complicated than the simple 
Newton-Lorentz equation, and they require computing 
b • and b x at each integration step. Still, they 
have the advantage that we can follow the overall mo- 
tion with relatively large time steps because we do not 



need to resolve the cyclotron motion. This reduces the 
cumulative error, as well as the total computation time. 

Figured] shows the solution of the guiding-center equa- 
tions for the same protons shown in Figure [2] under a 
dipolar magnetic field (Python source code provided in 
Supplement), 17 

It should be noted that the guiding-center equations 
are approximate because only terms first order in cy- 
clotron radius are used in their derivation. For parti- 
cles with larger cyclotron radii (higher kinetic energies), 
there may be a noticeable difference between guiding- 
center and full-particle trajectories (see Exercise[5]in Sec- 
tion EJ). 



V. THE SECOND AND THIRD ADIABATIC 
INVARIANTS 

The second adiabatic invariant is associated with the 
bounce motion, and it is calculated by integrating the 
canonical momentum over a path along the guiding field 
line: 



J2 = j> (jmv + qA) • ds, 



(24) 



where ds is the line element along the field line. The 
adiabatic integrals are evaluated in a "frozen" system: It 
is assumed that the drift is stopped, so the motion moves 
back and forth along a single guiding field line. 

Using Stokes' theorem, the second term can be con- 
verted to an integral over a surface bounded by the 
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FIG. 4. (Color online) Guiding-center trajectories for the par- 
ticles shown in Figure [2] The cyclotron motion is averaged 
out. 



bounce path 

j> qA-ds = q J V x A ■ da = q J B • da = (25) 

which is zero because the bounce motion goes along the 
same path in both parts of the cycle so that the enclosed 
area vanishes. Then, the second adiabatic invariant can 
be written as 



j 2 



rs m 2 

! / 7mv||ds, 



(26) 



where s is the path length along the field line, and s m i, 
s m 2 are locations of the mirror points where the particle 
comes back. 

At the mirror point the parallel speed v\\ vanishes so 
that v± = v. From the invariance of the magnetic mo- 
ment \i it follows that 
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(27) 




20 40 60 80 100 120 
time [s] 

FIG. 5. The second invariant values in time, calculated using 
the guiding-center trajectories in Fig. [H Lower curve is for the 
proton starting at 2 .Re distance, upper curve for the proton 
starting at 4 R e . 



where B m is the field strength at the mirror point. Sub- 



stituting v\ 



and solving for v\\ gives 



jmv\ 



B(s) 



ds = 2'ymvl. 



(28) 



The integral J2 is an adiabatic invariant in general (even 
if there are electric fields or slow time-dependent fields). 
If the speed is constant, I can be used as an adiabatic 
invariant. The integral / depends only on the magnetic 
field, not on the particle velocity, so it can be used to 
compute the drift path using the field geometry only (see 
exercises l8lfT0l in Section IVl]) . 

Figure [5] shows that the value of the second invariant /, 
evaluated using the guiding-center trajectories shown in 
Figured! stays constant in time. The integral is evaluated 
not using the definition of / in Eq. (|28j) . but using the 
dynamical form 



vtdt, 



7 ir 



(29) 



where the integral is evaluated over a half period. The 
limits of the integrals are determined by interpolation be- 
tween two points where the parallel speed changes sign. 
The values do not oscillate because the adiabatic invari- 
ant is calculated as an average over a cycle. 

The drift path is found by averaging the bounce mo- 
tion. In a dipolar field all drift paths are circular due 
to the symmetry of the field. The third invariant, associ- 
ated with the drift motion, is defined as an integral along 
the drift path: 



J3 = j> (jmv + qA) • dl, 



(30) 



where d£ is a line element on the drift path. 
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This can be written as 

J3 = <f> jmvddi ■ 



j) jmvddi + q J 



B da. 



(31) 



In the first term v&, the drift speed, is the magnitude 
of the expression given in Eq. (|2T|) . The second term is 
obtained by using Stokes' theorem as above. 

An order of-magnitude comparison shows that the first 
term of J3 can be neglected because it is much smaller 
than the second term: From Eq. (|2T|) the order of mag- 
nitude of the drift speed can be written as 

mv 2 B 

Vd ~2q^R> (32) 

where B is the typical field strength at the drift path and 
R is the typical distance from the origin. Similarly, from 
Eq. (|3j), the cyclotron radius has the order of magnitude 

mv 

P r o • 

qB 



(33) 



Then, the order-of-magnitude ratio of the terms in 
Eq. 43]} is 



(- 



(34) 



qBnR 2 q 2 B 2 R 2 

According to the adiabaticity condition Eq. (fl5| . p/R 
must be very small. Therefore the first term of Eq. (|3Tj) 
is ignored and we have 



^3 = q$, 



(35) 



where & is the magnetic flux through the drift path. The 
third adiabatic invariant is useful as a conservation law 
when the magnetosphere changes slowly, i.e., over longer 
time scales compared to the drift period. 

The use of three invariants gives more accurate results 
for the motion of particles over long periods. Numerical 
solution of equations of motion are less accurate because 
of accumulated numerical errors. Roederer 6 discusses 
in detail how drift shells can be constructed geometri- 
cally using the invariants (see Exercise fTOl in Sec. IVl]) . 
Furthermore, as three invariants uniquely specify a drift 
shell, the invariants themselves can be used as dynami- 
cal variables when investigating the diffusion of trapped 
particles. 20 



VI. EXERCISES AND ASSIGNMENTS 

This section lists some further programming exer- 
cises with varying difficulty. The code given in the 
supplement 17 can be modified to solve some of the ex- 
ercises. 

1. Uniform magnetic field. Follow charged parti- 
cles under a uniform magnetic field B = Bz where 
B = IT. Verify that the particles follow helices 
with cyclotron radius and frequency as given in 
Eqs. (j2j [3]). Experiment with particles with dif- 
ferent mass and charge values. 



2. Gradient drift. Consider a magnetic field given as 
B = (Ax + Bq)z. The field has a gradient in the in- 
direction, but no curvature. Set A = IT • m _1 and 
Bo = IT. Follow the trajectory of a particle with 
mass m = 1kg and q = 1C initialized with velocity 
v = lm • s _1 y at the origin. Note that the sideways 
drift arises from the fact that the cyclotron radius 
is smaller at stronger fields. 

3. Equatorial particles. Consider a particle in a 
dipolar magnetic field, located at the equatorial 
plane (z = 0) with zero parallel speed. As the field 
strength is minimum at the equator with respect 
to the field line, there is no parallel acceleration 
and the particle stays on the equatorial plane at all 
times. Using the dipole model, follow an equato- 
rial particle and verify that the center of the motion 
stays on a contour of constant £?, as implied by the 
conservation of the first adiabatic invariant. 

4. Explore the drift motion. Run the programs in 
the supplement to trace protons and electrons using 
the dipole model. Initialize particles with different 
energies, starting positions and pitch angles. Verify 
that electrons and protons drift in opposite direc- 
tions, and electrons have a much smaller cyclotron 
radius than protons with the same kinetic energy. 
Estimate the periods of bounce and drift motions 
and compare them with Eqs. ([7J [8]). 

5. Accuracy of the guiding-center approxima- 
tion. Simulate the full particle and guiding cen- 
ter trajectories with the same initial conditions and 
plot them together. Shift the initial position of the 
particle properly so that the guiding center runs 
through the middle of the helix. 

Repeat with protons with IkeV, lOkeV, lOOkeV 
and lMeV kinetic energies. At higher energies, the 
guiding-center trajectory lags behind the full par- 
ticle because the omitted high-order terms become 
more significant as the cyclotron radius increases. 

6. Different numerical methods. Solve the full 
particle and guiding center equations using differ- 
ent numerical schemes ^ 21 ! 22 such as Verlet, Euler- 
Cromer, Runge-Kutta and Bulirsch-Stoer. Verify 
the accuracy of the solution by checking the con- 
servation of kinetic energy and adiabatic invariants. 

7. Field line tracing. Plot the magnetic dipole field 
line starting at position (#o, Vo, zo). For any vector 
field u(r), a field line can be traced by solving the 
differential equation 



dr 



u 

R 



(36) 



where s is the arclength along the field line. 



8. Compute /. Compute the second invariant / 
(Eq. [28]) under a dipolar field for a guiding cen- 
ter starting at position (#o,2/o>0) and an equato- 
rial pitch angle a ec[ . The integral should be taken 
along a field line, which can be traced as described 
above. From the first adiabatic invariant one finds 
B m = B(xq, yo,0)/ sin 2 (a eq ) and the limits of the 
integral are found by solving B(s m ) = B m . 

9. Second invariant along the drift path. Pro- 
duce a guiding-center trajectory under the dipole 
field. By interpolation, determine the points 
(xi,yi,0) where the trajectory crosses the z = 
plane. Compute the second invariant / at these 
equatorial points and plot. Verify that the values 
are constant as shown in Fig. [5] 

10. Drift path tracing using the second invari- 
ant. Pick a starting location (#o?2/o>0) and mir- 
ror field B m , and evaluate the second invariant 
I(xo, yo, B m ) as described above. Compute the gra- 
dient V/ = d x I£ + dyly numerically using central 
differences: 

d x I ~ [H x o + S, 2/ 0j B m ) - I( x o — ^ 2/o, B m )] (37) 
d y I « [I(xo, yo + S, B m ) - I(x , y - 5, B m )] ,(38) 

where S is a small number (e.g. 0.01R e ). 

The second invariant is constant along the drift 
path, so for a finite step (Ax, Ay), it holds that 
d x IAx + d y IAy = 0. Use this relation to trace suc- 
cessive steps along the drift shell. This method is 
more accurate than following a particle or a guiding 
center. 

11. The double-dipole models The dipole ceases 
to be a good approximation for the magnetic field 
of the Earth as we go farther in space. The double- 
dipole model, although unrealistic, introduces a 
day-night asymmetry that vaguely mimics the de- 
formation of the magnetosphere by the solar wind. 
It can be used to capture some basic features of 
particle dynamics in the outer magnetosphere, if 
only qualitatively. 

The model has one dipole (Earth) at the origin, 
pointing in the negative z-direction, and an image 
dipole at x = 20R e . If both dipoles are identical, 
the magnetic field is given by: 

B(x, y, z) = B dip (x, y, z) + B dip (x - 20R e , y, z) (39) 

where B d i p (x, y, z) is given by Eq. (j6j). 

The domains of each dipole are separated by the 
plane x — 10 R e . This plane simulates the magne- 
topause, the boundary between the magnetosphere 
and the solar wind. For slightly better realism, the 
image dipole can be multiplied by a factor larger 




y[R e ] 



FIG. 6. A proton with 200keV kinetic energy, initialized 
on the right edge at (7R e ,7R e ,0) with 60° pitch angle in a 
double-dipole field. 

than 1 so that the magnetopause becomes curved. 
Also, the two dipoles can be tilted by equal and 
opposite angles with respect to the Sun-Earth line 
(x-axis), to simulate the fact that the dipole mo- 
ment of the Earth is tilted. 

(a) Starting at various latitudes, plot the mag- 
netic field lines of Eq. ([39]) on the x — z plane. 
Observe the compression of field lines on the 
day side and extension on the night side. Note 
that no field line crosses the x = 10R e plane. 
Multiply the image dipole term by 1.5 and re- 
peat. 

(b) Follow several guiding-center trajectories 
starting position x between —7R e and — 10i? e ? 
y = z = 0, and pitch angles between 30° 
and 60° (smaller pitch angle creates a longer 
bounce motion). With small pitch angles, the 
particle should come closer to Earth on the 
day side. Repeat with a pitch angle of 90°. 
Now the particle goes away from the Earth on 
the day side. Explain these observations using 
the conservation of first and second adiabatic 
invariants. 

(c) The double-dipole field can break the second 
invariant for some trajectories. The reason 
of this breaking is that the field strength has 
a local maximum on the dayside around the 
equatorial plane. Particles with sufficiently 
small mirror fields are diverted to one side of 
the equatorial plane because they cannot over- 
come this field maximum, as seen in Fig. [6] 
Start an electron guiding center at position 
xq = — lOi^e, yo = zq = with kinetic energy 
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K = lMeV with an equatorial pitch angle 70° 
and follow its guiding center for 1000 seconds 
with time step 0.01. On the dayside the tra- 
jectory will temporarily move above or below 
the equatorial plane. Using the method used 
in Fig. [5j plot the second invariant / versus 
time. The second invariant will be constant 
between breaking points, but its value will dif- 
fer from the initial value. The reason is that 
near the breaking points bounce motion slows 
down and the adiabaticity condition does not 
hold. However, the first invariant is not bro- 
ken. 

This phenomenon, named drift-shell 
bifurcation^ can be one of the causes 
of particle diffusion in the magnetosphere. 




12. Magnetotail current sheet. On the tail region of 
the magnetosphere, magnetic field lines are heavily 
stretched, and a sheet of current is flowing through 
them. 24 The field in the magnetotail can be repre- 
sented by the simple form: 

B = £ ^x + £ n z if |z|<l 

= + 5 n z otherwise (40) 

In this problem we set Bq = 10, B n = 1 and 
d = 0.2. The field lines trace parabolas on the 
x — z plane, which can be seen by integrating the 
equation dx/B x = dz/B z . The parameter d is the 
scale of the current sheet thickness. The truncation 
of the field at z = 1 simulates the finite size of the 
tail region. 

The field vector points to opposite directions on 
both sides of the equatorial plane z = 0. When 
a charged particle is released from above, it moves 
toward the weaker region near \z\ =0 where the 
adiabaticity condition does not hold. The helix be- 
comes a "serpentine orbit" that moves in and out 
of the equatorial plane. The chaotic dynamics of 
these orbits is extensively studied . 25 ! 26 

Figure [7] shows the three types of orbits that can ex- 
ist in such a mode h 26 i 27 "Speiser orbits" approach 
the equatorial plane and later go beyond z = 1 and 
leave the tail region. "Cucumber orbits" alternate 
between helical and serpentine orbits. These do not 
form closed orbits because of the breaking of the 
first invariant at the equatorial plane. "Ring or- 
bits" alternate between oppositely-directed fields; 
they do not have a helical section. 

(a) Evaluate for \z\ < 1 and determine the 
direction of gradient-curvature drift. 

(b) By trial and error, find initial conditions that 
create the types of orbits shown in Fig. [71 




FIG. 7. (Color online) Types of orbits created by a particle 
with mass m = 5 and charge q = 1 near a current sheet, (a) 
Speiser orbits of transient particles, (b) Cucumber orbits of 
quasitrapped particles and (c) Ring orbits of trapped parti- 
cles. Note the different scales of axes. 

VII. CONCLUDING REMARKS 

Space plasmas provide many case studies which, after 
proper simplification, can be used in the undergraduate 
physics curriculum. We have presented one such case, the 
basic theory of charged-particle motion under the dipole. 

This paper focuses on visualization and concrete com- 
putation, with the hope that students will modify or 
rewrite the code to run their own numerical experiments 
on particle motion in magnetic fields. In my opinion, nu- 
merical simulations provide at least two important ped- 
agogical benefits: First, even if the required analytical 
tools are beyond the students' level, they can use sim- 
ulations to obtain a qualitative understanding. Second, 
the process of coding the simulation forces students to 
understand the problem at a basic and operational level. 

The main body of this article or the exercises can be 
incorporated in lectures, or they can be given as ad- 
vanced assignments to interested students. A natural 
place for this subject is a course on electromagnetism 
and/or plasma physics. When the basics are introduced, 
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the instructor can discuss related subjects such as plasma 
confinement, radiation belts, or space weather. 

In advanced mechanics courses, adiabatic invari- 
ants are usually presented with an abstract formalism. 
Charged particle motion provides a natural and concrete 
case where adiabatic invariants are relevant and indis- 
pensable. 

The subject can also be incorporated in courses on 
computational physics. Accuracy and stability of differ- 
ent numerical integration schemes may be presented us- 



ing charged particle motion. The widely separated time 
scales of the motion would be a challenge for most of the 
schemes. 
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